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Abstract 

A  numerical  technique  capable  of  simulating  blade-scale  compression  system  flow  instabilities 
over  times  scales  spanning  tens  of  rotor  revolutions  is  presented.  Simulations  of  stall  inception, 
growth  to  fully-developed  rotating  stall,  and  evidence  for  hysteresis,  secondary  instabilities,  and 
other  nonlinear  phenomena  are  presented.  Signal  processing  techniques  for  flow  asymmetry 
characterization  are  discussed  in  the  context  of  obtaining  low-order  representations  of  the  flow 
disturbances  with  the  ultimate  goal  of  active  stall  suppression. 

Keywords:  compressor  stall,  instabilities,  bifurcations,  computational  fluid  dynamics,  co¬ 
herent  structures 


1  Introduction 

Models  for  understanding  the  aerodynamics  of  compression  system  stall  inception  and  growth 
of  these  disturbances  to  fully-developed  rotating  stall  have  existed  for  over  forty  years.  A 
feature  common  to  many  of  these  models  is  that  the  power  input  to  the  flow  field  is  modeled 
by  some  forcing  function  to  the  equations  of  fluid  motion  [1,  2].  This  “indirect”  approach  is 
favored  over  direct  numerical  simulation  because  of  the  wide  time  scale  range  of  the  dynamic 
phenomena  (several  revolutions  for  stall,  10-100’s  for  a  surge  cycle)  and  the  complexity  of  direct 
computational  fluid  dynamic  simulations  of  multistage  compression  systems. 

Recently,  numerical  techniques  which  represent  the  force  imparted  to  fluid  flow  fields  as  forc¬ 
ing  functions  to  the  Navier-Stokes  equations  [3]  (replacing  boundary  conditions  in  the  flow  and 
pressure  fields)  have  been  developed  for  simulating  systems  with  fluid/flexible  structure  inter¬ 
actions  [4] .  Because  the  speed  with  which  fluid  dynamic  simulations  can  be  performed  increases 
with  decreasing  complexity  of  the  domain  over  which  the  computations  are  performed,  these 
techniques  can  also  improve  the  efficiency  of  these  computations.  Exploitation  of  the  interplay 
of  this  numerical  technique  and  the  communication  patterns  of  massively-parallel  processors  is 
a  key  element  in  making  the  multistage  compressor  simulations  possible  in  this  work. 
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The  modeling  approach  presented  in  this  study  is  unique  in  the  field  of  turbocompressor 
simulation  in  that  not  only  do  the  forcing  functions  depend  on  spatial  position,  but  they  are 
also  computed  solely  from  the  equations  of  fluid  motion  and  the  geometry  of  the  compressor 
blading.  An  important  aspect  of  this  work  is  studying  the  predictive  capabilities  of  such  a 
model  in  simulating  the  inception  of  rotating  stall  as  the  stall  margin  is  crossed,  and  following 
the  growth  of  these  instabilities  to  fully-developed  rotating  stall.  Simulations  of  rotor  blade  wake 
momentum  defect  growth  to  a  fully  separated  flow  field  blocking  a  blade  passage  will  be  shown 
to  be  consistent  with  the  conceptual  model  of  compressor  stall,  originally  proposed  by  Iura  and 
Rannie  [5].  Stall  cell  propagation  rates  of  approximately  1/2  rotor  speed  were  observed  in  our 
simulations.  The  stall  cell  propagation  mechanism  put  forth  by  the  cited  study1  corroborates 
with  our  simulations. 

The  simulation  results  presented  in  this  paper  focus  on  a  two  stage  compression  system 
model.  Each  rotor  row  consists  of  11  blades  and  each  stator  has  14  blades.  The  number  of  blades 
was  chosen  to  exceed  the  minimum  number  of  rotor  blades  involved  in  a  propagating  modal  stall 
cell:  Longley  [6]  estimated  that  the  minimum  wavelength  of  a  modal  stall  disturbance  would 
span  eight  blades,  and  this  theoretical  prediction  is  supported  by  experimental  studies  with 
a  seven  rotor  blade  machine  [7]  which  resulted  in  stall  events  consistent  with  the  continuum 
models  of  Moore  and  Greitzer  [1],  Naturally,  smaller-scale  events  should  then  be  captured  in 
our  design. 

Understanding  the  spatial  structure  of  the  stall  cells  is  important  for  interpreting  signals 
obtained  from  a  finite  number  of  measurement  probes  in  experimental  systems  and  is  particularly 
important  for  designing  stall  controllers  based  on  direct  suppression  of  these  flow  disturbances. 
Early  studies  of  the  spatial  structure  of  stall  cells  assumed  single  sin-shaped  modes,  and  studies 
which  built  on  this  initial  work  decomposed  the  stall  cells  into  modes  of  a  truncated  Fourier 
series.  We  will  present  preliminary  results  on  obtaining  a  optimal  (empirically  determined)  set  of 
basis  functions  from  our  detailed  simulations,  a  technique  which  has  the  potential  of  generating 
extremely  low  order  simulators  which  retain  virtually  all  of  the  fidelity  of  the  original  simulations. 


2  The  Numerical  Simulation  Technique 

The  compression  system  model  developed  in  this  paper  is  based  on  2D,  incompressible  flow 
through  an  annulus,  and  so  we  begin  by  writing  the  Navier-Stokes  equations  in  conservative 
form  and  separating  the  right-hand-side  into  terms  which  are  linear  and  nonlinear  in  the  velocity 
field  v  to  obtain 


O  -j 

-tt  =  -V  •  (vv)  +  — V2v  -  Vp  +  f  =  N(v)  +  L(v)  -  Vp  +  f  (1) 

at  Re 

over  the  unit  square.  What  follows  is  an  outline  of  the  numerical  technique  used  to  integrate 
(1)  over  one  time  step  At.  The  time  step  of  all  simulations  performed  in  this  study  was  fixed  at 
0.001  rotor  revolutions.  A  staggered  arrangement  of  velocity  component  discretization  points 
was  used  in  the  finite-difference  calculation  of  the  spatial  derivatives  [8].  Pressure  and  forcing 
function  scalar  field  values  were  defined  at  the  cell  centers  using  a  500  x  500  evenly  spaced  array 
of  cells.  Dimensionless  rotor  speed  is  1.0  and  the  Reynolds  number  based  on  the  circumferential 
length  was  set  at  Re  =  10, 000. 

At  the  start  of  each  time  step,  given  a  velocity  field  vn  at  time  tn  and  following  the  operator 
splitting  method  of  Karniadakis  and  co-workers  [9],  we  predict  an  interim  velocity  field  v  from 

1The  mean  flow  field  diverted  around  the  blockage  results  in  blades  to  one  side  experiencing  increased  angles  of 
attack,  thus  inducing  flow  separation,  and  blades  to  the  other  side  experiencing  the  opposite  effect,  resulting  in  overall 
stall  cell  propagation  in  the  former  direction 
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< —  Rotation 

Figure  1:  Compression  system  geometry  on  the  discretization  grid  points.  Mean  axial  flow  is  from 
bottom  to  top.  Rotors:  1st  and  3rd  blade  rows;  stators  are  2nd  and  4th.  Vy„,  marks  the  axial 
velocity  measurement  locations. 

a  fourth-order  Runge-Kutta  (explicit)  step: 

V  —  vn  1 

=  g  (k0  +  2kl  +  2k2  +  k3) 

with 

k0  =  N(v”)  k2  =  N(v”  +  kj  At/2) 

kx  =  N(v"  +k0Al/2)  k3  =  N(v"  +  k2  At). 

Adding  the  contributions  of  the  pressure  and  forcing  fields  gives 

=  -Vp  +  fb  +  ft  with  ft  •  ey  =  (2) 

where  acts  only  in  the  direction  ey  of  the  mean  axial  flow  and  is  applied  over  the  entire  flow 
field  as  a  substitute  for  the  atmospheric-to-plenum  pressure  rise.  V  is  the  mean  axial  outlet 
velocity  component,  and  7  is  the  throttle  position  in  the  orifice  equation. 

The  spatially  resolved  forcing  function  fb  is  computed  to  impart  the  equivalent  amount  of 
force  that  would  otherwise  be  “felt”  by  the  fluid  if  the  fluid/blade  boundary  condition  actually 
existed.  If  x  is  the  vector  of  k  discretization  points  defining  an  individual  rotor  or  stator  blade, 
v(x,;)  is  the  interim  velocity  field  interpolated  to  the  discretization  point  x,:,  and  v*  is  the  blade 
velocity, 

k 

ST  W  V*  -  v(Xj) 

2^<l(x,;,xi)fb(xj)  =  - — - . 

j= 1 

Solving  this  system  of  linear  equations  for  the  individual  forces  fb  is  numerically  efficient,  since 
the  array  of  <5  values  only  must  be  inverted  once.  This  numerical  simplification  is  possible  since 
the  arrangement  of  forcing  points  on  each  blade  remains  constant,  and  because  the  blades  are  far 
enough  from  each  other  so  that  their  force  fields  do  not  overlap.  The  numerical  distance  function 
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Figure  2:  Stall  inception  after  throttle  position  change  from  7  =  0.3  to  7  =  0.2. 


S(xi,xj)  used  above,  in  calculating  v(x*),  and  for  spreading  fb  to  the  pressure  discretization 
point  locations  is  discussed  in  Peskin  [4] . 

Continuity  is  enforced  by  taking  the  divergence  of  (2)  and  solving  the  resulting  Poisson’s 
equation  for  the  static  pressure  field  p : 


V2p 


V  •  v 
At 


+  Vfb. 


(3) 


Solving  the  (spatially)  discretized  form  of  (3)  amounts  to  solving  250,000  coupled,  linear,  alge¬ 
braic  equations  -  a  large  number  resulting  from  the  fine  mesh  size.  This  means  iterative  solution 
methods  must  be  used  over  a  direct  method,  such  as  Gaussian  elimination.  Because  we  have 
replaced  what  would  otherwise  be  a  complicated,  time-dependent  geometry  with  forcing  func¬ 
tions  defined  over  a  simple  domain,  we  can  take  advantage  of  relatively  non-memory  intensive 
and  naturally  parallelizable  techniques,  such  as  the  Conjugant  Gradient  Method.  A  simple  do¬ 
main  means  less  use  of  routed  and  more  nearest-neighbor  communications  in  massively-parallel 
processors,  and  so  we  have  found  significant  advantages  to  performing  our  simulations,  written 
in  CM  Fortran,  on  a  Thinking  Machines  Corporation  CM5. 

The  final  step  in  computing  the  total  update  also  benefits  from  fast  linear  equation  solvers 
since  we  chose  to  use  an  implicit  integrator  for  calculating  the  contribution  of  the  linear  portion 

old): 


v  ”  1  -  v 


At 


L(v”+1). 


This  concludes  the  discussion  of  the  numerical  integration  technique. 
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3  Simulation  Results 

One  of  the  hallmarks  of  compression  system  aerodynamic  instabilities  is  the  hysteresis  associ¬ 
ated  with  compressor  stall.  In  a  noise-free  situation,  we  find  that  as  the  compressor  throttle  is 
closed,  the  plenum-to-atmospheric  pressure  rises  while  the  mean  axial  flow  decreases  until  the 
axisymmetric  flow  field  becomes  unstable.  After  this  point  is  crossed,  small  perturbations  grow 
rapidly  into  large- amplitude  rotating  stall  cells.  Reopening  the  throttle  does  not  necessarily 
reverse  this  transition:  the  multistability  brought  about  by  the  nonlinear  nature  of  the  instabil¬ 
ities  can  lead  to  nonrecoverable  stall  and  operating  points  which  are  unstable  with  respect  to 
finite  amplitude  disturbances.  Noise  clouds  the  stall  inception  picture,  since  small  amplitude 
inlet  flow  field  disturbances  can  be  amplified  by  the  impending  instability  and  may  give  rise  to 
the  frequently  cited,  long-lived,  small  amplitude  stall  precursors. 

Preliminary  simulations  were  performed  to  find  the  approximate  location  of  the  stall  incep¬ 
tion  point.  The  compression  system  simulator  was  then  brought  to  equilibrium  at  a  throttle 
opening  corresponding  to  7  =  0.3,  a  “uniform-flow,”  locally  asymptotically  stable  operating 
point  (we  place  the  term  uniform-flow  in  quotes  since  the  velocity  field,  even  in  the  inlet  and 
exits  ducts,  is  never  truly  uniform  due  to  flow  around  the  compressor  blading).  The  throttle 
was  then  closed  through  the  range  where  the  flow  was  found  to  become  unstable  to  a  final  value 
of  7  =  0.2.  The  time  trace  of  Fig.  2  shows  the  onset  of  stall  soon  afterwards.  Note  that  the 
disturbance  grows  to  a  large-amplitude  rotating  stall  cell  after  only  a  few  rotor  revolutions, 
growing  only  out  of  the  natural  flow  asymmetry  induced  by  the  rotor-stator  blade  interactions. 
The  early  state  of  this  flow  disturbance  is  illustrated  in  Fig.  5. 

After  a  period  of  disturbance  growth  to  fully  developed  rotating  stall,  the  throttle  was  re¬ 
opened  to  7  =  0.3  (recall  that  this  throttle  position  corresponds  to  a  locally  asymptotically 
stable  uniform  flow  operating  point).  After  reaching  equilibrium,  we  found  that  the  rotating 
stall  persisted,  a  clear  indicator  of  multistability  (hysteresis).  See  Fig.  3  for  a  close-up  view  of 
the  stalled  flow  field.  Furthermore,  plotting  the  mean  exit  pressure  versus  time  (Fig.  4)  shows 
a  secondary,  time-dependent  oscillation.  Because  the  frequency  is  so  much  lower  than  the  fre¬ 
quency  with  which  stall  cells  interact  with  individual  rotor  blades,  we  take  this  as  evidence  of 
secondary  bifurcations  taking  place  along  the  stalled-flow  equilibria  locus  born  off  the  axisym¬ 
metric  flow  solution  branch.  Because  of  the  incompressible  formulation  of  this  simulation,  the 
exit  duct  does  not  act  as  a  mass  storage  region  and  so  this  is  not  a  surge  bifucation,  but  results 
from  the  birth  of  a  modulated  traveling  wave.  See  Adomaitis  and  Abed[10]  for  more  details  on 
the  bifurcation  behavior  of  these  systems. 

4  Signal  Processing 

One  of  the  primary  motivations  for  this  research  was  to  develop  a  compression  system  flow 
instability  simulator  with  no  built-in  bias  towards  the  structure  of  the  stall  inception  and  fully 
developed  flow  field  instabilities.  This  approach  would  give  a  model  capable  of  resolving  both 
modal  and  localized  disturbances,  and  so  would  provide  an  effective  tool  for  developing  numerical 
techniques  for  characterizing  the  flow  asymmetry  and  testing  signal  processing  techniques  for 
identifying  such  flow  features  in  an  experimental  rig.  To  this  end,  consider  the  problem  of 
finding  a  set  of  spatially-dependent  trial  functions  </>,:(:c)  from  which  a  time-dependent  linear 
combination  can  be  computed  to  represent  the  the  axial  component  of  flow  field  ( Vy ): 

Vy(x,t )  =  a0(t)tl>o(x)  +  (nU)ih  (x)  +  •  •  •  +  am(t)^m(x). 

Taking  =  1  means  00  represents  the  mean  axial  flow  I4iean,  and  all  1/7  for  i  =  1,2,... 
represent  spatially-varying  modes,  so 

v(x,t)  =  ai{t)tpi(x)  +  . . .  +  am{t)ipm{x).  (4) 
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Figure  3:  Stable  flow  field  relative  to  rotor  1 ,  illustrating  stalled  and  unstalled  blade  passages  (to 
the  right  and  left,  respectively)  for  7  =  0.3.  The  stall  cell  propagates  to  the  right  in  this  reference 
frame.  Note  the  region  of  reversed  flow  upstream  of  the  stall  cell. 
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Figure  4:  Plenum-to-atmospheric  pressure  trace  during  rotating  stall  for  7  =  0.3.  The  6-rev 
oscillations  (a  much  lower  frequency  than  that  generated  by  stall  cells  passing  from  blade- to-blade) 
indicates  secondary  instabilities. 


Ultimately,  we  will  consider  low-order  approximations  to  the  entire  2D  flow  field,  but  this 
study  will  be  limited  to  a  ID  “slice”  across  the  circumferential  coordinate  at  the  axial  position 
of  the  velocity  measurement  probes  (see  Fig.  5).  Taking  snapshots  of  the  instantaneous  axial 
velocity  component  of  the  flow  field  deviations  at  this  location  facilitates  comparisons  to  previous 
studies,  such  as  Moore  and  Greitzer[l]  where  the  flow  field  disturbance  just  upstream  of  the  inlet 
guide  vanes  was  modeled. 

Because  finite  difference  approximations  were  used  to  discretize  (1),  each  snapshot  taken 
during  the  stall  inception  transient  of  Fig.  3  consists  of  a  long  (500  element)  vector  of  velocity 
values.  Each  of  these  vectors  corresponds  to  points  in  time  separated  by  approximately  1/3  rotor 
rev  intervals,  giving  a  total  of  M  =  43  snapshots  collected  during  the  12  rotor  rev  transient. 
Prior  to  applying  the  numerical  technique  to  be  discussed,  we  estimate  the  phase  angle  of  the 
flow  field  disturbances  and  adjust  all  of  snapshots  to  synchronize  the  phase  angles  so  that  the 
minimum  velocity  point  all  coincide  at  the  Vy§  probe  position. 

The  time  evolution  of  the  flow  field  is  governed  by  set  a  well-behaved  differential  equations 
and  so  there  is  some  degree  of  correlation  between  the  snapshots.  This  means  the  snapshots  are 
not  necessarily  linearly  independent,  and  so  an  efficient  method  for  representing  the  snapshot 
vectors  is  to  first  determine  an  orthogonal  basis  spanning  the  snapshots  and  then  decompose  the 
snapshots  using  the  basis  vectors  as  discretized  trial  functions.  One  can  use  the  Gram-Schmidt 
orthogonalization  procedure  to  obtain  the  basis  vectors,  but  an  optimal  basis  is  determined  by 
a  computationally  efficient  implementation  of  the  Proper  Orthogonal  Decomposition  method 
discussed  in  Sirovich  [11]. 

If  the  vector  Vi  represents  the  deviation  of  the  axial  flow  component  Vy  from  the  mean  axial 
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Figure  5:  Flow  field  axial  velocity  component  deviation  in  the  inlet  duct  at  the  axial  position  of 
the  velocity  measurement  probes  Vyn. 


flow  Vmean  of  snapshot  i,  we  can  construct  the  array  of  inner  products: 


(vi,Vi) 

(v2,vi) 

(VM,V  l) 

and  so  determine  the  optimal  set  of  eigenmodes  spanning  the  original  set  of  snapshots  using  the 
eigenvectors  Et  of  (5).  Each  element  Ejyi  gives  the  contribution  of  snapshot  j  to  eigenmode  i, 
and  so  the  empirical  eigenfunctions  are  computed  by 

M 

i’dx)  =  ^E^v^x). 

j= i 


(r?i  j  Vo ) 
(V-2,V2) 

(VM,V2) 


(vi,vM) 

(v2  ,vm) 

( Vm,Vm ) 


(5) 


These  eigenfunctions  are  shown  in  Fig.  6. 

With  each  eigenmode  ipi  is  associated  an  eigenvalue  A,:  which  corresponds  to  the  probability 
of  finding  that  eigenfunction  in  the  time-dependent  flow  field  from  which  the  snapshots  were 
extracted: 

A,;  =  [(Ui  +  (V2  ,^i)2  +  •  •  •  + 


An  optimal  set  of  trial  functions  should  capture  as  much  energy  as  possible  in  as  few  modes  as 
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0.25 


Figure  7:  Contribution  of  the  first  four  eigenmodes  to  the  snapshots  collected  during  stall  inception. 


possible,  and  this  holds  true  for  the  example  presented  since 

'  Ai  =  0.7179  ' 

A2  =  0.1212 
A3  =  0.0679 
A4  =  0.0313 
A5  =  0.0189 
A6  =  0.0131 
A7  =  0.0095 
A8  =  0.0052 


showing  that  nearly  96%  of  the  energy  is  captured  in  the  first  five  inodes. 

These  empirically  determined  eigenmodes  can  now  be  used  to  decompose  snapshots  of  the 
dynamically  changing  flow  field  to  quantify  growth  rates  of  these  modes  as  the  stall  cell  develops. 
Projecting  snapshots  of  the  flow  field  onto  the  basis  functions  gives  the  values  of  the  mode 
amplitude  coefficients  of  (4);  results  are  shown  in  Fig.  7. 

One  can  reverse  the  decomposition  procedure  and  reconstruct  the  flow  field  from  mode 
amplitudes  and  the  empirical  eigenfunctions,  such  as  shown  in  Fig.  8.  This  is  more  than  simple 
data,  compression,  since  we  see  that  the  first  mode  may  provide  a  more  efficient  better  way 
of  representing  the  flow  asymmetry  than  the  first  term  in  a  trigonometric,  wavelet,  or  other 
theoretical  basis  function  series  since  it  captures  72%  of  the  disturbance  energy.  Analogous  to 
the  Galerkin  procedure  used  by  [1],  we  can  use  these  empirical  eigenfunctions  as  a  basis  for 
discretizing  the  original  equations  of  fluid  motion  giving  nonlinear  mode-amplitude  equations 
which  should  predict  the  bifurcations  and  dynamics  in  the  neighborhood  of  where  the  snapshots 
were  collected.  In  fact,  the  eigenmodes  of  Fig.  6  can  be  used  without  modification  for  discretizing 
the  model  of  Moore  and  Greitzer  [1], 
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Circumferential  Coordinate 


Figure  8:  Axial  flow  deviation  from  mean  (left)  and  flow  deviation  reconstructed  from  the  first 
eight  modes  (right)  corresponding  to  the  flow  field  at  rotor  revs=6  of  Fig.  2. 


The  flow  field  distubances  can  be  faithfully  reconstructed  for  more  developed  rotating  stall 
cells  and  the  small  amplitude  disturbances  at  stall  inception.  Figure  9  is  a  stall  inception 
disturbance  (c.f.  Fig.  5,  the  true  disturbance)  reconstructed  from  the  same  eight  eigenmodes  of 
Fig.  6,  showing  the  initial  flow  disturbance  is  localize  to  several  blades.  It  is  interesting  to  see 
that  the  blade  passage-width  oscillations  grow  in  amplitude  to  the  side  of  the  main  flow  deviation 
corresponding  to  the  higher  angle  of  attack  (the  circumferential  coordinate  range  spanning  0.8 
to  1.0  and  then  continuing  from  0.0  to  0.2)  and  are  suppressed  on  the  other  side  (0.4  to  0.8),  in 
accordance  to  the  stall  cell  propagation  mechanisms  discussed. 
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